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1 Introduction 

Spherically-symmetric models are of interest in the study of the nonlinear inhomogeneities 
of the universe for two reasons. First, spherical symmetry is often a reasonable working 
approximation for the inhomogeneities in the real universe. Second, spherical symmetry 
allows to exactly solve the Einstein's equations, a highly nontrivial task. Recent research 
has focused mainly on the Lemaitre model [1] (see [2-5] for recent contributions) which 
describes the dynamics of a spherically-symmetric perfect fluid, and great attention has 
received its pressureless limit, usually named the Lemaitre- Tolman-Bondi (LTB) model [1, 
6, 7], which can also include a nonvanishing cosmological constant. LTB metrics have been 
used to describe the large-scale inhomogeneities of the late universe as, for example, in Swiss- 
cheese [8-15], void [16-36] and inhomogeneous [37-39] models for apparent acceleration (see 
[40, 41] for recent reviews). Besides being useful in understanding the role of voids in the 



- 1 - 



process of structure formation [42, 43], the LTB model is also suitable to study the dynamics 
of spherical collapse in an expanding universe, black- hole formation included [44-48]. 

In the present paper we extend the Lemaitre model to the case of n decoupled and 
non-comoving perfect fluids with general equations of state. Our main result is the set of 
3 + 2n exact equations governing the dynamics of the model, which we write in terms of 
physically meaningful quantities and express in a concise and transparent way. Many are 
the possible applications, both at early and late times. At late time, one can study the 
evolution of overdensities and underdensities in a universe where dark energy is not the 
cosmological constant and, in particular, can be inhomogeneous (see, for example, [49, 50] 
and references therein). This possibility may also induce an inhomogeneous variation of 
fundamental constants, as for instance the fine structure constant if a coupling between the 
dark energy and the electromagnetic field is allowed (see, for example, [51-55]). Moreover, 
dark matter and baryons can be described as two separate fluids, the latter possibly featuring 
pressure. As stressed in [5], the introduction of non-dust equations of state can have a non- 
negligible effect on the cosmological models and this should be taken into consideration 
while interpreting cosmological datasets. This is, however, a non-trivial task due to the 
increased number of free parameters, and such difficulties increase when one considers n 
fluid components with the consequence that more astronomical and astrophysical input is 
needed to constrain the available parameter space. Finally, at early times the contribution of 
radiation can be included, which may be relevant both for the understanding of the evolution 
of the standard post-inflation inhomogeneities and for the correct modeling of the very large 
underdensities typical of void models [34] . 

This paper is organized as follows. In Section 2 we will go through all the details of our 
formalism, and in Sections 3 and 4 we will present the numerical results for the case of two 
non-comoving dust components in a flat ACDM universe and for spherical collapse in the 
presence of dark energy with negligible speed of sound, respectively. Then in Section 5 we will 
compare our findings to previous work dealing with exact solutions. Finally we will give our 
conclusions in Section 6. In Appendix A we discuss the expansion tensor and in Appendix 
B the case of the Lemaitre metric, possibly in a non-comoving frame. The definitions of the 
functions used to model the numerical example of Section 3 are given in Appendix C. 

2 The model 

2.1 Metric and Einstein tensor 

A spherically-symmetric metric can be written as 

y'2 

d s 2 = - e 2X dt 2 + — dr 2 + Y 2 dQ 2 , (2.1) 

1 + 2E 

where the lapse function e A , the scale function Y and the curvature function E depend on 
the coordinate time t and coordinate radius r, d$7 2 = d9 2 + sin 2 6 d(j) 2 and we have set c = 1. 
A prime denotes partial derivation with respect to the coordinate radius r, whereas a dot 
denotes partial derivation with respect to the coordinate time t. Comma and semicolon 
signs will not denote derivatives. We will use a reference frame comoving with the arbitrary 
four-velocity field u^, which has components = (e _A , 0, 0, 0). Moreover, we assume 
that 7'/0 and E > —1/2 in order to always have a regular g rr > 0. We will discuss more 
general evolutions and shell crossing in forthcoming work. 
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The nontrivial components of the Einstein tensor G a p = R a p — \Rg a p for the met- 
ric (2.1) are then: 



G tt = e 2X 



G tr = 2e x X'H A - E, 



Y' 

diy 



Y' 2 ( o , 1 + 2B . 2£ 

(E 1 1 -I- 2E 
AR + A A + H R H A -yjy- -YJ^XF 

where we have defined radial and angular expansion rates: 

_ VTT2E ^TOT _ e~ x Y' 
Hr ~—i ~ ~Y 1 ~ ~ dl ' (2 - 2) 

and also radial and angular acceleration rates: 

A R = ^P^^fP = ^ " 2^1^ + El - E d2 - e~ x XH R , (2.4) 

1 d 2 Y e~ 2X Y 
4 

respectively, and also the following auxiliary quantities: 



= ^^ = ^-e" A A^, (2.5) 



e- x E e~ 2X E Y' Y" E' X" 

Edl ~lT2E> Ed2 ~\T2E' • F -Y~Y 7+ lT2E +X+ Y- 

Note that for A = and E = the usual LTB expressions are recovered. 

In the previous equations d/dr T f = uf { d a = e~ x d/dt is the derivative with respect to 
the proper time of the comoving observer. As we show in Appendix A, the components of the 
acceleration of u" { are a r>I f = X' and a t>T t = 0. We see therefore that, for a geodesic reference- 
frame velocity field, the lapse function A depends only upon time and can be rescaled such 
that A = and d/dr r { = d/dt. From Eqs. (2.2-2.3) it follows then that in a non-geodesic 
reference frame the comoving observers measure with different proper times. Later we will 
identify the reference- frame velocity with the velocity field of one of the fluid components, 
and it will turn out that the acceleration a r r f = A' is sourced by pressure gradients which 
push the observers out of the freely-falling geodesies. 

2.2 Einstein's equations and energy-momentum tensor 

From the Einstein's equations G a $ = n T a p we can obtain four independent dynamical equa- 
tions. We choose to form two of them with the following combinations: 

YVGu-Y^G* = {e . 2XYf2 2EY) , = Y--Y'T lt -Y>YT tr ^ 

r I U rr — 1 1 Utr _ / -2A-iat>2 r, p V V 1 I l rr — X Y ±tr 



YY 2 - 2EY) = k—^jj—^ , (2.7) 



y ,2 /(l + 2£) V > -Y' 2 /{l + 2E) 
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where k = 8nG and T a @ is the energy-momentum tensor for an ideal fluid source comprised 
of n decoupled components: 

n 

T a/3 = Y^ T ^, (2.8) 

i=l 

where the energy- momentum tensor of the i:th perfect fluid component is: 

Tf = Pl u^+ Pl hf, (2.9) 

where uf , p\ and pi are the four- velocity field, energy density and pressure, respectively, 
and hj = g a/3 + ufu^ is the projection tensor on the hypersurface orthogonal to uf. The 
isotropic pressure is related to the energy density by p = w p, where the equation of state 
parameter w = w(t,r) is assumed to be a general function of t and r (see, for example, [56] 
for the case of a static fluid with anisotropic pressure). In the chosen coordinate system it is 

uf = 7*(e-\ Vi, e , 0, 0), (2.10) 

where vi c is the coordinate comoving peculiar velocity of the i:th component relative to the 
reference frame. The proper peculiar velocities and the gamma factors are instead given by 

y/2 1 
v Ip = YT2E V ^ c 7 ^ = 1 - v 2 ' (2 ' U) 

respectively. The covariant components of the total energy-momentum tensor are then: 

n 

T tt = e 2X Pi + w iH ~ w i] » ( 2 - 12 ) 
i=l 

Ttr = -^^e^piil + wjvi,^?, (2.13) 

^ i=l 

Y,Pi [{I + Wi)vl pl f + W t ] , (2.14) 



1 + 2E . 



T e0 = T^/ sin 2 6 = Y 2 Y,Pi- (2-15) 



i=i 



As the third dynamical equation we will consider the combination coming from — G\ + 
G r r -\- G d d + G^, which gives the generalization of the acceleration equation: 

VIT2E (Y' \ kA ; . . 

A = A R + 2A A -^— Orff— +.FJ =--Y J Pi{l + Zwi)i (2.16) 

^ ' 1=1 

where we used the fact that the acceleration scalar of the reference-frame velocity is a r f = 
^ X y 2E A' (see Appendix A). The term proportional to a r f gives a "spurious" contribution 
to the acceleration and vanishes if we use a geodesic reference frame, in which the total 
acceleration is Ar + 2Aa, similarly to LTB models. 

Finally, the last independent equation will be simply the Gt r = & Tt T component, which 
reads: 

e- x E k Y . . o e~ x Y 

y+2~e = 2vmE^ p * {1+m)VM + vmE ° rf • (2 - 17) 

z=i 
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Eq. (2.17) shows that the evolution of the curvature is due to two distinct causes. The first is 
the effect of having an inhomogeneous multicomponent fluid and goes to zero in the FLRW 
limit where the peculiar velocities vanish or if only one fluid is present and its reference 
frame is adopted. 1 It is sourced by the energy flux in the radial direction: we remind indeed 
that the curvature function E may be interpreted as the total energy of a given shell (see 
Eq. (2.28)). The second contribution to E is due to the fact that, generally, we are using a 
non-geodesic reference frame (a r f 7^ 0) in which the total energy of a shell r is not conserved. 

2.3 Conservation equations 

As we are considering an ideal fluid source comprised of n decoupled components, each i:th 
component will be conserved individually, i.e., V = 0, where V Q denotes the covariant 
derivative. For each i:th component we obtain the two following equations 2 where we omit 
the index i and multiply for clarity's sake by e x and v c iX_ 2 e > respectively: 



v c Y' 2 
1 + 2E 



e x V a T at = {p + p)j (6 + v p a) + v cl 2 {p + p)' + e~ x pv 2 pl 2 + e" W = 0, (2.18) 

ar /„ 1 /„,2a 1 „, \ 1 „.2_,2„-A/„ , „v , „/„, „,2 , „, „/„,2„,2 



V Q T ar = {p + p) 7 (v 2 p @ + v p a) + v 2 pl 2 e-\p + p) + p'v cl 2 + v c p'j 2 v 2 = 0, (2.19) 



where the expansion scalar and the acceleration scalar a have been calculated in Appendix 
A. By taking the combinations e x V a T at - v c ^ E V a T ar and e x V a T at - ^-V a T ar it is 
possible to obtain the relativistic energy-conservation and Euler equations, respectively, for 
non-comoving fluids: 

^ = - @ (p + p), (2.20) 

^ = -a(p + p), (2.21) 

where, analogously to the total derivative with respect to the proper time r, we defined the 
convective derivative with respect to a along the four-acceleration a^ 1 : 



d „„ _ x d d 



dr _„ % = 7e -»_ +7 „ c _, (2 ^ 2) 



Note that in the rest frame of the i:th fluid component the derivative with respect to r 
reduces to {—gtt)^ l ^ 2 d/dt and the derivative with respect to a to g r }^ 2 d jdr. Eq. (2.21) is 
actually a contracted form of the standard Euler equation, which reads W v d v p = —a^ (p+p). 
From Eq. (2.21) it is clear that the four-velocity of a dust (p = 0) component is geodesic. 

By identifying the reference- frame velocity uf f with the i:th fluid-component velocity 
u? we can obtain the rest-frame expressions of Eqs. (2.20-2.21) by simply setting v c = 0: 

e- A Pi = -e rf ( ffi + pj), (2.24) 
Pi = -a r ,rt(pi + Pi) , (2-25) 



1 Note, however, that the cosmological constant never sources E. 

2 The other two equations relative to the angular components identically state the spherical symmetry of 
the energy-momentum tensor. 



- 5 - 



respectively, where a r r f = A' and a t r f = (see Appendix A). The last equation shows that 
the pressure gradients are the cause for the non-geodesity of the reference frame relative to 
the i:th fluid component 

Finally, because of the (twice contracted) Bianchi identity it holds V^C" = and so: 

n 

5^V M 3f = 0. (2.26) 

i=i 

Consequently, two conservation equations will not be independent and can be discarded. 
Equivalently, we may replace, if convenient, any other two dynamical equations with the 
latter two extra conservation equations. 

2.4 Gravitating mass 

The effective gravitating total mass F is given by 

2GF = H A Y 3 - 2E Y , (2.27) 

where G is the gravitational constant. This quantity corresponds to the mass in the Schwarzschild's 
metric if the metric of Eq. (2.1) was matched to the latter at a given radius r. See Refs. [4, 
57, 58] for more details. It is interesting to note that Eq. (2.27) may be rewritten as: 

1 f dY\ 2 GF . 

from which follows that the curvature function E can be interpreted as the total energy per 
unit of mass of a shell, and that to it contribute the kinetic energy per unit of mass and the 
potential energy per unit of mass due to the total gravitating mass up to that shell. Note that 
thanks to spherical symmetry one is able to define a potential energy also in cases far away 
from nearly Newtonian ones and that the potential energy is related to the curvature [7]. 

In giving the initial conditions it will be useful to divide the total mass F into the 
separate masses of the i:th fluid components: 



i=l 



Note, however, that this division is irrelevant for the dynamics as only the total F enters the 
evolution equations. By substituting Eqs. (2.12-2.14) and (2.27) into the (combinations of 
the) Einstein's Eqs. (2.6-2.7) one obtains the equations for F and F' , which using Eq. (2.29) 
read: 

n n 

e- X ^F i = -4vrY 3 £ Pi 7 f [H A (^ p + w { ) + S A (1 + Wi )] , (2.30) 

i=l i=l 
n n 2 

• 1 • 1 C J 

1=1 1=1 ' 

where the angular "spatial" expansion rate S A = v c y- has been introduced in Appendix 
A. Note that the Einstein's equations give Eqs. (2.30-2.31); therefore to write Eq. (2.29) 
we have put the constants of integration to zero. Note also that the evolution of the i:th 
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mass Fi is sourced by the gravitational influence of pressures and velocities of the other fluid 
components, differently from the case of the conservation equations (2.20-2.21) which apply 
to the fluid components individually: decoupled fluids are indeed still gravitationally coupled. 
In the case of only one component Eqs. (2.30-2.31) have a simple interpretation as discussed 
in Appendix B. 

2.5 Full set of equations 

The total number of independent equations is generally 4 + 2n — 2, where the first term is the 
number of nontrivial Einstein's equations, the second comes from the conservation equations 
and the last from the Bianchi identity. By including the definition of mass in Eq. (2.27) we 
have a total of 3 + 2n equations. We choose to use all the 2n conservation equations, two of 
which will replace the acceleration equation (2.16) and the evolution equation (2.30) for the 
total F. We now list for clarity the full set of 3 + 2n exact equations governing the dynamics 
of a spherically symmetric inhomogeneous model with a multi-component ideal fluid: 

tfi = ^E^ + ^, (2.32) 

i=l 

e- x E k Y 2 e~ x Y 

y+2~e = 27rrH^ M1 + ^ Kp74 + 7rrH arf ' (2 ' 33) 

n n 2 

Y, F i= 47ry3 X> — \ h a0- + w *) v Ip + Sa(1 + vl p Wi)] , (2.34) 



i=l i=l 
dp 



dn 
dpi 
do; 



®i(Pi+Pi) i = l,...,n, (2.35) 

ai(pi+pi) i = l,...,n, (2.36) 



where the unknown functions are Y, E, the total F and n copies of pi and t> j jC . The reference- 
frame velocity can be taken to be geodesic so that a rt = A = 0. Alternatively, one can 
identify the reference- frame velocity with the velocity of the i:th fluid component, = u?, 
with the consequence that vj c = and A becomes an unknown function in its stead. The 
equations of state Wi have to be given as an external input. From Eq. (2.32) it seems clear that 
the dynamics is similar to the FLRW one, the complication being an array of interconnected 
equations defining the evolution of gravitating mass and curvature. 

From Eq. (2.32) one may obtain the age of the universe at a given coordinate radius r. 
For example, in the case of expansion (positive root, Y > 0) it is: 



t A (r,t) = t - t B (r) = [ dte~ X Y 

Jt B (r) 



2GF(r,t) , 



-1/2 

(2.37) 



which can be computed once the dynamical equations (2.32-2.36) are solved. The quantity 
ts is the so-called bang function also present in LTB models. 

2.6 Initial and boundary conditions 

Initial conditions. We will give initial conditions at the time t < to, where to is the present 
time. First we fix the residual gauge freedom in the coordinate r by setting Y(r, t) = a(t)r, 
where a(t) is the scale factor of the external FLRW model to which we will match the 
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inhomogeneous metric (see below). If one is not interested in embedding the metric, a(t) 
may be taken as an arbitrary number. Then, E(r,t), F{r,t) and n copies of pi{r,t) and 
Vi,c(r,t) are needed. 

The curvature function E may be specified by demanding a homogeneous age of the 
universe t& at t, which gives an additional constraint between the total F and E. From 
Eq. (2.37), in the approximation that F(r,t) ~ F(r,F) and E(r,t) — E(r,t) for t < i , we can 
obtain: 

rY(r,t) 

t A (r,i) = i-t B (r)~ / dye- x ^ 
Jo 

which can be used without solving the dynamical equations to demand at t the homogeneous 
big bang, ie(r) = 0. 

Note then that in the set of Eqs. (2.32-2.36) there is not a dynamical equation for F, 
which can always be found by integrating Eq. (2.34) for a constant t once the boundary 
condition F(0,t) = is given. In particular, Eq. (2.34) relates the initial conditions for F, 
Pi and Vi^ c , which cannot be given independently, as it is intuitively clear. It is important to 
point out that while the initial conditions can be given individually for the different fluids 
by considering Eq. (2.34) as n independent equations, this will not be true at later times 
for which only the total F matters. In other words, one cannot obtain the individual Fi by 
integrating the corresponding i-piece of Eq. (2.34) and then combining them to obtain F. 
This also means that for t > i there is not a meaningful way to define the individual Fi, 
which are not physically relevant. We can instead define the invariant mass M as shown, in 
the example of Section 3, by Eq. (3.3). 

If we identified the reference- frame velocity with the velocity of the i:th fluid component, 
u£f = uf, we may have to give initial conditions for A (and not rr i c which is identically zero). 
Note, however, that A only appears in the acceleration equation which we have discarded. 
As with F , it is therefore enough to give a boundary condition (see below) in order to obtain 
with Eq. (2.25) the lapse function at i for any r. 

We have now all the functions and their r-derivatives on the hypersurface of constant 
i. The next step is to evolve the model forwards in time along the worldlines of constant r, 
which is done by simultaneously integrating the first-order in t differential equations listed in 
(2.32-2.36). While the t-derivatives for Y, E and pi are evident, Vi c enters only through the 
expansion and acceleration scalars. Finally, A(r, t) and F(r, t) are found by integration for 
the needed constant t. The equations that we have discarded thanks to the Bianchi identity 
may be used to cross-check the accuracy of the numerical integration. 

Boundary conditions. In discussing the necessary boundary conditions we will partic- 
ularize the analysis to a metric which is matched to an external FLRW model at some 
comoving coordinate radius r&. We will give boundary conditions at the initial time: the 
evolution equations will automatically maintain them at later times. Similarly to the embed- 
ding of an LTB model, the curvature and the integrated gravitating masses have to match 
the FLRW values, E(rt,,i) = —\kr\ and Fi(rb,t) = F° nt (rt,,t), respectively (in the chosen 
gauge the scale function is already matched as Y(r,t) = a(t)r). The latter condition for 
the case of Vi >c (r,t) = explicitly demands that /9° ut (t) = ■% /q 6 pi{f ,t) f 2 df , as it is easy 

b 

to see from Eq. (2.34). According to the Darmois-Israel boundary conditions [59] the lo- 
cal density may be discontinuous through a timelike surface (constant r&), but the pressure 
must be continuos Pi(rb,t) = p° ut (i). Furthermore the peculiar velocities have to vanish, 
Vi, c (rb,t) = 0, as also the lapse function, \(rb,t) = 0, as we want to describe the FLRW 



2GF{r,t) 



+ 2E(r, t) 



(2.38) 
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model in a geodesic reference frame where proper and coordinate times coincide. Even if not 
necessary it may be desirable to match smoothly the background, i.e., pi(r{,,t) = p° ut (t) and 
Pi(n,t)=p' i (r b ,t) = X'(r b ,t)=vl c (r b ,t)=0. 

Finally, we discuss the regularity conditions at the center of the inhomogeneous metric. 
Because we do not want a singularity at the center it is -Fj(0, t) = E(0, t) = 0. In particular, 
near the origin it is F$ oc r 3 and E oc r 2 , as one can see from Eq. (2.27) where Far 
and a finite Ha 7^ is assumed (see also [37]). The latter also implies Y oc r and so 
Y(0,t) = Y(0,t) = 0. Also vanishing at the center must be the velocities, Vi tC (0,t) = 0, and 
the pressure gradient, 2^(0, t) = 0, which implies A'(0, t) = 0. It seems again desirable to 
have densities and velocities smooth at the origin, p^(0, t) = v[ c (0, t) = 0. 

2.7 Light propagation 

The metric (2.1) implies that for radial null worldlines it is 

dt dr , Y' , 
— = ±^e~ x . 2.39 
du du \J\ + 2E 

where u is an affine parameter. Let us consider a photon emitted at t(u) in a coordinate 
time v(u) (emitting source and observer are taken to be in the rest frame). Expanding the 
right-hand side of Eq. (2.39) around t(u) and keeping terms up to first order in v{u) we get 

dt dv dr x Y' dr Y' f TT x -\ .„ 

^ + ^ = -3^ e ^^^=-j — , [H R -e A A u, 2.40 

du du du y 7 ! + 2E du y/l + 2E V / 



which, using Eq. (2.39) for dt/du, gives 

'H R -e- A \\. (2.41) 



1 dv dr Y' 1 , , _ A 



v du du \J\ + 2E 

The redshift is defined as 3 



e x °u - e x v 



(2.42) 



where the subscript refers to values evaluated at the affine-parameter value uq at the 
observer's position. The derivative of z with respect to u is then 

dz sdr ( dX 1 du\ ,„ . dr ( Y' TT \ .„ 

— = -(1 + z )— — + -— = - 1 + z)— X' - H R . 2.43) 

du K 'du\du uduj v ' du \ y/l + 2E J 

Using (2.43) and remembering the expression for the rest-frame acceleration scalar (see Ap- 
pendix A), Eq. (2.39) gives a pair of equations 

dt _ 1 e~ x dr _ 1 1 y/1 + 2E 

dz 1 + z a r f — Hr ' dz 1 + z a r f — i^R Y"' 

which determines the coordinates t and r on the light cone as functions of the redshift, and 
can be put in the following gauge invariant form: 

dr rf = _L 1 ; ±E = L_ 1 ; (2.45) 

dz 1 + z a r f — ///J, ' 1 + z arf — ' 

3 This definition is different from the one in [5]. The authors of [5] agree, however, that (2.42) is correct. 
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where r p is the proper radial distance. In the ALTB model where E = a r f = 0, the previous 
equations reduce to the familiar form 

£ = Y —, *= (2.46) 

^ (i + «)y <fe (i + z)y 

3 Two non-comoving dust components in a flat ACDM universe 

As a first application of the formalism introduced in this paper, we will consider an inhomo- 
geneous sphere embedded into a flat ACDM universe. The inhomogeneities will be given by 
two non-comoving dust components, one of which may be interpreted as baryonic matter and 
the other as dark matter. The reader should keep in mind, however, that these results are 
merely illustrative of the multi-fluid model here presented. This solution becomes the usual 
ALTB model if the density of the second dust component is set to zero or if their relative 
velocity is set to zero. 

We will use the reference frame comoving with the component labelled by "M" , so that 
the relative Euler equation gives A = 0. The other conservation equation reads: 

™=-{H R + 2H A ) = -^f, (3.1) 
Pm JrJa 



and can be directly integrated to give: 

p M (r,t) = J R (r,t)J A (r,t) = Y 2 (r,t)Y'(r,t) y/l + 2E{r,t) 
p M {r,t) J R {r,t)J A {r,t) ^/l + 2E(r,t) Y 2 (r,t)Y'{r,t) ' 



(3.2) 



where i is the initial time. We will label the other dust component with "TV" and use the 
simpler notation 7^r = 7, vn, p = v p , vn,c = v c as the corresponding M quantities are trivial. 
About A note that the corresponding conservation equations (2.35-2.36) trivially state that 
PA is uniform and constant. 

3.1 Initial and boundary conditions 

First we set the flat ACDM model by fixing the present-day expansion rate to Hq = lOO/i 
km s _1 Mpc -1 with h = 0.7 and the matter density parameter to flatter = 0.3. We then 
solve the background equations to find the evolution of the scale factor a(t) and so the initial 
time t corresponding to z = c ^jy — 1, which we will fix illustratively to z = 5. By redshift z 
we will always mean the redshift relative to an observer in the background model, which will 
differ, for example, from the redshift relative to the observer at the center of the spherical 
inhomogeneity. The inhomogeneous sphere will have a comoving radius of r& = lOO/i -1 Mpc. 

Next we give initial and boundary conditions at t. For the scale function we set Y(r,t) = 
a(t)r. Then we have to give initial conditions for F(r, t) = Fm{t) + F^{r) + F\(r). For A it 
is simply F\{r) = 4^a 3 (t)r 3 /OA- We stress again that the division of the total F into separate 
components is meaningful only at the initial time as only the total gravitating mass F enters 
the evolution equations. Note, however, that it is possible to define individual invariant mass 
components Mj. For instance, for the M component it is: 

AttY 2 Y' 
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from which it is easy to verify using Eq. (3.1) that it is correctly Mm = (while F^O). 

The remaining background matter density is split between the M and TV" components 
according to a fraction q such that it is pffiit) = #Pmatter(*) an d PjvH*) = (1 — q) Pmatter(*)- 
We choose the profile of the M component to be: 

out /TN = 1 + Oi{r, r tM ,d M i,M2, o>mi,mv I 3 - 4 ) 

Pm W 

where the 5 n function is given in Appendix C and 5 m = 5 mi + <$M2 is the contrast in matter 
density between the center of the spherical inhomogeneity and the background value of p ^ ■ 
This inhomogeneous profile satisfies the smoothness conditions p' M (0,t) = p' M (rb,t) = and 
continuously matches the background density Pm = Pm*(^)- Moreover, the matching 

requires that Fm{t) has the corresponding background value and so the following condition 
has to be satisfied: 

PM = P PM(f,i)r 2 df, (3.5) 

from which follows that it has to be 5 mi 5m2 < 0. 5 mi > will give an overdensity, while 
5 mi < an underdensity. For a given model we will fix the quantities 5mi,M2 and olm\,M2 
and we will get analytically the radius rt M by demanding Eq. (3.5). The gravitating mass is 
then simply Fm{t) = 47ra 3 (f) f£ pM(f,t)r 2 df, which has also an analytic expression. 
Next we need jFjv(r), which we will take to be: 

F N (r) = F^\r) +5 3 (r,r tN ,5 NhN2 ,a NhN2 ) , (3.6) 

with 5n2 = —8m so that Fjy(r) is correctly matched to the background value of F^ ut (r) = 

We can now find the initial curvature function E by demanding that the universe has the 
same age t for any r, i.e., from Eq. (2.38): tA_(r, t) = t. The approximation from neglecting 
the time dependence of F and E for t < i is very good as the cosmological constant is 
negligible at early times. 

Next, we have to set the peculiar velocities for the N component with respect to the M 
component. As there is not pressure, the case v c (r, t) = gives trivially the dynamics of the 
ALTB model with a density profile that is the average of the M and N density profiles. We 
choose the initial velocity profile to be: 

v c (r,t) = 5i(r,r tv ,5 v i jV2 ,a v i, v2 ) (3.7) 

with 5 v2 = —5 v x so that there are no peculiar velocities at the center and at the border. The 
contrast 5 v2 will give the maximum peculiar velocity in c units. Finally we can compute the 
initial density profile for the iV component, which from Eq. (2.34) is: 

PN(r,t)= _ _ 2/ M ) rT , (3.8) 

where Ha is given by Eq. (2.32) evaluated at the initial time. This terminates the necessary 
initial conditions. 

3.2 Evolution 

We will now show the exact general relativistic evolution of the system presented in the 
previous Section for the case of a central underdensity. 4 The left panel of Fig. 1 shows the 



The actual parameters we use are q — 0.7, ftjui = om2 = 0.4, 8mi = —0.3, Sm2 = 0.05, r tjv = 2r(,/3, 
ajvi = 0.1, Qjv2 = 0, S N1 = -S N2 = 0.1F£ ut (r b /2), r tv = 2r b /3, a vl = a v2 = 0, S v2 = -5 vl = 0.5 • 10" 2 /a(t). 



-11- 



rbkg=Y/a (Mpc/h) f _ Fa 




20 40 60 80 100 20 40 60 80 100 



r (Mpc/h) r bk g (Mpc/h) 

Figure 1. Left panel: evolution of the scale function Y for the times corresponding to the background 
rcdshifts of z =5 (red), 3, 1 and (blue). In particular, this plot shows the relation between the 
coordinate radius r and the background comoving radius rbkg = Y/a which will be used in all the 
other plots. Right panel: Evolution of the gravitating mass F for z =5 (red), 3, 1 and (blue) with 
the contribution due to A subtracted and normalized to unity at hole radius r&. The time evolution is 
purely an effect of the peculiar velocities between the two dust components. See Section 3.2 for more 
details. 



evolution of the scale function Y: in the interior region the scale function is expanding faster 
because of the less matter and consequently negative curvature present, as also shown by 
the spatial Ricci scalar 1Z in the left panel of Fig. 2 (see Appendix B for its definition). All 
the other figures will be plotted with respect to the background comoving radius defined as 
r bkg = Y(r,t) /a(t). This is a good definition as the curvature function is E <C 1 (see right 
panel in Fig. 2) and thus the proper radial distance is r p = f J[_^2E ~ ^ ' ' 

The evolution of the scale function Y is governed by the gravitating mass F and the 
curvature function E, as shown by Eq. (2.32). The evolution of the gravitating mass F is 
plotted in the right panel of Fig. 1. We have subtracted the time-dependent component due 
to A, which is -Fa = Y 3 PA, m order to show the genuine time dependence of the gravitating 
mass due to the dust components. As there is no pressure, the evolution displayed in the plot 
is due to the peculiar velocities between the two components. In particular the gravitating 
mass at half the radius is decreasing in accordance with the outgoing flow shown in the right 
panel of Fig. 4. 

The evolution of the curvature function E, is plotted in the right panel of Fig. 2. This 
is again a genuine effect of the peculiar velocities between the two dust components, as a 
cosmological constant never sources the evolution of the curvature function E. We remind 
that E can be interpreted as the total energy per unit of mass of a given shell (see Eq. (2.28)), 
and one can see from the plot that energy seems to flow outwards, again in agreement with 
the peculiar velocity field shown in the right panel of Fig. 4. As said before, the left panel 
of Fig. 2 shows the Ricci scalar 1Z, which seems to follow a profile which is an average of the 
density contrasts of the two dust components, displayed in Fig. 3. We also remind that the 
curvature function is related to the Euclidean average of 1Z, as pointed out in Appendix B. 

Both the density contrasts of the two dust components show a realistic evolution, as one 
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Figure 2. Left panel: evolution of the spatial Ricci scalar 1Z plotted as the dimensionless combination 
lZa 2 r% for z =5 (red), 3, 1 and (blue). Its profile is closely related to the density profiles of the 
two dust components displayed in Fig. 3. Right panel: evolution of the curvature/energy function E 
for z =5 (red), 3, 1 and (blue). The plot shows an outgoing flow of energy in accordance with the 
peculiar velocity field shown in the right panel of Fig. 4. See Section 3.2 for more details. 
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Figure 3. Evolution of the M component local density (left panel) and N component local density 
(right panel) for z =5 (red), 3, 1 and (blue). Both the density contrasts show a realistic evolution, 
with overdense regions contracting and becoming thin shells (mimicking structures), and underdense 
regions becoming larger and deeper (mimicking voids). Note in particular how the N component 
induces gravitationally a peak in the density of the M component. See Section 3.2 for more details. 



can see from Fig. 3. Overdense regions start contracting and they become thin shells (mim- 
icking structures), while underdense regions become larger and deeper (mimicking voids), and 
eventually they occupy most of the volume. An interesting feature of this multi-component 
model is that the gravitational interaction between the two fluids is evident. The profile of 
the M component (left panel) develops indeed a peak exactly where the N component peaks: 
that is, an overdensity of one components drags in the collapse the other component. This 
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Figure 4. Left panel: evolution of the peculiar velocities Vm = Y{Ha — H out ) of the M component 
with respect to the background for z =5 (red), 3, 1 and (blue). Velocities are increasing as the 
underdensity becomes emptier. Right panel: evolution of the peculiar velocities of the N component 
with respect to the M component for z =5 (red), 3, 1 and (blue). The gravitational attraction 
between the two components reduces their relative velocities. See Section 3.2 for more details. 
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Figure 5. Same as Fig. 4 but with the velocity of the peak subtracted. These plots show that 
the velocity flow near the peaks in density is natural: matter is falling towards the overdense regions 
which indeed grow with time as illustrated in Fig. 3. See Section 3.2 for more details. 



can be checked by decreasing the matter density allocated to the N component, with the 
effect that the the peak at 90/i _1 Mpc in the M component disappears. 

Finally, in Fig. 4 we show the peculiar velocities of the two components. The left panel 
shows the peculiar velocities of the M component with respect to the background, which 
are given by vm = Y(Ha — H out ). This definition again relies on the fact that r p ~ Y. 
As expected these peculiar velocities are growing with time as the underdensity becomes 
emptier and emptier (see also the left panel of Fig. 1). Opposite is, however, the case of 
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the peculiar velocities of the N component with respect to the M component (right panel of 
Fig. 4), which decrease with time: the gravitational attraction between the two components 
dumps their relative velocities. The plots of Fig. 4 are with respect to the "reference frame" 
of the background and of the center of the underdensity, as the velocities go to zero at r = 
and r = r&. In Fig. 5 we re-plot the same curves, but with the velocity of the peak in 
density subtracted. Fig. 5 therefore better shows the velocity field close to the overdensity, 
which appears to be natural, as matter is falling towards the peak in the density from both 
directions. 

Finally, this model maintains the scale invariance valid for the simpler LTB model [9] if 
the initial peculiar velocities are properly scaled: v c (r,t)/ri, = v c (r, t)/^ where tilde marks 
the quantities relative to an inhomogeneous patch of different radius. 



4 Collapse in the presence of dark energy with negligible speed of sound 

As a second application of our general model we will consider an inhomogeneous sphere 
embedded into a flat wCDM universe. The inhomogeneities will be given by a pressureless 
matter component and by a dark energy source with constant equation of state w ou t < —1/3. 
This solution becomes the usual ALTB model if w = w Q ut = — 1 and the Lemaitre model 
(possibly in non comoving coordinates) if the density of the dust component is set to zero. 

We will use the reference frame comoving with the matter component, which we label 
with M, so that A = and Eqs. (3.1-3.2) hold. We will label the dark energy component 
with u w" and use the simpler notation w w = w, p w = p, j w = 7, v WiP = v p , v W;C = v c 
as the corresponding M quantities are trivial. As we are discussing a fluid with negative 
pressure we cannot consider adiabatic perturbations for which w(r,t) = w ou t: this would 
indeed give a negative speed of sound {c 2 s = w out ) which would cause an unstable growth of 
perturbations [60]. We will consider instead the following r-dependent non-adiabatic equation 
of state: 



w(r,t) = w ont 
which implies the following speed of sound: 



Pw{r, t) 



(4.1) 



c 2 s = ^ = aw(r,t), (4.2) 

so that for a < the speed of sound is positive (the adiabatic speed of sound is recovered 
with a = 1, see e.g. Ref. [61]). 

Particularly interesting is the case a = 0, for which the pressure gradients are absent 
and the speed of sound vanishes, with the consequence that matter and dark energy collapse 
together following the geodesic flow. The possibility of a dark energy fluid with negligible 
speed of sound has been investigated in the literature under various assumptions, and gener- 
ally requires a non canonical scalar field like /c-essence, as the standard quintessence models 
with canonical scalar fields always have c s = 1 (see, for example, [49, 50, 62-64] and refer- 
ences therein). Apart from theoretical considerations, the fact that the speed of sound of dark 
energy may vanish opens up new observational consequences. Indeed, the absence of dark 
energy pressure gradients allows instabilities to develop on all scales, also on the small scales 
where dark matter perturbations become non-linear. We expect, therefore, dark energy to 
modify in a detectable manner the growth history of dark matter halos not only through its 
different background evolution but also by actively participating to the structure formation 
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process, in the linear and non- linear regime [65-68]. The observational consequences of a 
clustering dark energy have been indeed extensively studied. See, for example, [69-76] for 
the impact on large-scale structures and cosmic microwave background. 

Finally, we would like to point out that, as we will use the matter reference frame, 
the speed of sound given in Eq. (4.2) will not be in general in the dark energy rest frame, 
which is usually used in order to have a gauge independent expression for the speed of sound. 
However, as we will not introduce initial peculiar velocities, if c s = then the matter and 
dark energy rest frames will always coincide because of the absence of pressure gradients. 
The latter will also be approximately true for the case of nonzero but negligible speed of 
sound as the peculiar velocities developed by the dark energy fluid will be very small (see 
the right panel of Fig. 9). 

4.1 Initial and boundary conditions 

We fix the background model by setting h = 0.7 and fiji/o = 0.3 as in Section 3. About the 
equation of state we will study the following three cases: uwt = —0.8 and a = in Section 
4.2, w out = —0.8 and a = — 10~ 10 in Section 4.3 and w ou t = —1.2 and a = — 10~ 10 in Section 
4.4. The inhomogeneous sphere will have a comoving radius of r& = lO/i -1 Mpc and we will 
give initial conditions at the time i corresponding to the background redshift z = 1000. We 
set again the scale function at initial time to Y(r, t) = a(t)r. 

We have to give initial conditions for F(r,t) = Fm{t) + F w (r). We model the initial 
matter inhomogeneity as: 

f m(t) / t\ 1 - tanh((r - r - r offset )/2Ar) 

F^r)- l + dM[t) 1 + tanh(r /2Ar) ' (43) 

where Fffi = ^a 3 (t)r 3 p^ t (t) is the background gravitating mass, 5m is the density contrast, 
and the parameters tq and Ar characterize respectively size and steepness of the density 
profile. The parameter r g- set is used to translate the profile in order to have a smooth origin. 
We will adopt the following values: tq = r fj se t = 0.25 n,, Ar = 0.3 vq and <$jtf(i) = 1.5 • 10~ 3 , 
that is, we have a central overdensity. The initial matter density is then: 

Next, we can find the initial curvature function E by demanding that the universe has 
the same age i for any r, i.e., from Eq. (2.38): t^ir, t) = t. As F w {r) Fm(t) at t, this 
condition will be independent of the dark energy initial conditions. In particular at the 
initial time the model is very close to a (dust) LTB model, so that it is possible to find E 
analytically by means of the following expression valid for 5m "C 1 [77]: 

E(r,t) * -l(a(i)Ho» t (t)r) 2 - l) . (4.5) 

These initial conditions are free from decaying modes [78]. 

Finally, we take the dark energy component to have no initial peculiar velocities, which 
means that it will have, with respect to the background, the peculiar velocities YAH of the 
matter component. Therefore, in order to have initial conditions without decaying modes in 
the dark energy component, we have to set the initial profile as: 

F w (r) 1 + Sw{i) 1 - tanh((r - g , - r oS t )/2Ar) 



F° ut (r) w 1 +tanh(r /2Ar) 
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Figure 6. Left panel: evolution of the scale function Y for the times corresponding to the background 
rcdshifts of z = 1000 (red), 100, 10, 2, 1, and (blue). In particular, this plot shows the relation 
between the coordinate radius r and the background comoving radius rbkg = Y/a which will be used 
in the other plots. Right panel: evolution of the peculiar velocities vm = Y(Ha ~ H out ) of the M 
component with respect to the background. Velocities are directed inwards and increase in magnitude 
with time as the central overdensity becomes denser. See Section 4.2 for more details. 



where F° ut = ^ L a 3 (f)r 3 / o° ut (t) and initial contrast is: 

Ut)= ! + ™° Ut M*), (4-7) 

1 - 6W ou t 

which is valid during matter domination for sub- Hubble perturbations and c^Cl [79]. The 
initial dark energy density is then: 

^ = 45§F- < 48 > 

With this we have all the necessary initial conditions. 
4.2 Evolution for zero speed of sound 

The left panel of Fig. 6 shows the evolution of the scale function Y for w ou t = —0.8 and a = 0: 
in the interior region the scale function is expanding slower because of the more matter and 
consequently positive curvature present. All the other figures will be plotted with respect 
to the background comoving radius defined as rbkg = Y(r,t)/a(t), which as explained in 
Section 3.2 is a good definition of distance. The right panel shows the peculiar velocities of 
the M component with respect to the background. As there are no initial decaying modes, 
the peculiar velocities are natural, with the matter falling towards the peak in the density, 
as shown by the left panel of Fig. 6. Moreover, the velocities grow with time as the central 
overdensity becomes denser and denser. 

In Fig. 7 the evolution of the density contrasts of the matter and dark energy components 
is shown. The evolution, as remarked above, is realistic. In particular the dark energy 
component closely follows the evolution of the matter component, as one should expect from 
a fluid with vanishing speed of sound. Because of the absence of pressure gradients peculiar 
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Figure 7. Evolution of the matter contrast (left panel) and dark-energy contrast (right panel) for 
z = 1000 (red), 100, 10, 2, 1, and (blue). Both the density contrasts show a realistic evolution, 
with the overdense region becoming denser and denser. Note, in particular, how the dark energy 
component follows the evolution of the matter component, with the difference of a slower growth rate 
due to the nonzero pressure. See Section 4.2 for more details. 



velocities between the two components never develop (if initially zero) and the two fluids 
both follow the same geodesic flow. 5 As it is clear from the plot, however, the growth rates 
are different. This has to be expected because the negative pressure generally slows down the 
evolution of perturbations. To be quantitative, we have plotted in right panel of Fig. 8 the 
evolution of the ratio of the contrasts, S w /Sm, which we normalized to the value valid during 
matter domination, pjg^ 1 ■ As one can see, the agreement is perfect when ~ 1 (see 
left panel of Fig. 8), showing that the chosen initial conditions are indeed free from decaying 
modes. During the nonlinear dark-energy-dominated phase the ratio departs as expected 
from the matter-domination prediction, and in particular becomes different in overdense and 
underdense regions, as illustrated in the right panel of Fig. 8 where the blue line is relative 
to rbkg = l/i _1 Mpc (overdensity) and the green line to rbkg = 6/i _1 Mpc (underdensity) . 
This could be due to the fact that even though the two components have similar profiles, in 
the overdense region the matter component reaches the nonlinear evolution earlier (the dark 
energy is actually almost linear at the present time). 

It is also interesting to compare the evolution to the case of no dark energy. We find 
that the matter perturbation grows to a maximum contrast of about ~ 90. Alternatively, we 
find that an initial perturbation of 5m (t) = 1.1 • 10 -3 grows to the same present-day contrast 
as the initial perturbation of <5a/(*) = 1-5 ■ 10~ 3 studied in this Section. 

Finally, let us comment on the fact that because the speed of sound is zero, there are 
no characteristic length scales associated to the dark energy clustering and the spherical 
collapse remains independent of the size of the object [65]. The characteristic length scale 
associated to the dark energy clustering is indeed the sound horizon scale, L s = a J c s dt/a, 
which vanishes for c s = so that clustering takes place on all scales. In other words, this 
model maintains the scale invariance valid for the simpler LTB model [9]. 

5 For the same reason the curvature is time independent, as it is easy to see from Eq. (2.33). 
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Figure 8. Left panel: redshift evolution of the matter and dark energy density parameters. Right 
panel: evolution of the ratio of the contrasts, 5 w /Sm, normalized to the value valid during matter 
domination, l*^ ut ■ The blue (dot-dashed) line is relative to rbkg = lft> _1 Mpc and the green (solid) 
line to rbkg = 6/i Mpc, see Fig. 7. See Section 4.2 for more details. 




r bkg (Mpc/h) r bkg (Mpc/h) 

Figure 9. Left panel: evolution of the speed of sound c 2 s for the times corresponding to the 
background redshifts of z = 1000 (red), 100, 10, 2, 1, and (blue). Right panel: evolution of the 
peculiar velocities of the dark-energy component with respect to the matter component. See Section 
4.3 for more details. 



4.3 Evolution for c s ~ 10 -5 

We will now consider the model of the previous Section (io ut = —0.8) with a = — 10 -10 , 
which corresponds to a speed of sound of order c s ~ 1CT 5 , as one can see from the left panel 
of Fig. 9. This latter value implies a sound horizon of order L s ~ 80/i _1 kpc, much less than 
the inhomogeneity scale considered in this example, so that the overall dynamical evolution 
shown in Figs. 6-8 is basically unchanged. The new features of this case are pressure gradients 
and peculiar velocities. In order to understand their behavior it is useful to work out the 
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Figure 10. Evolution of the matter contrast (left panel) and dark-energy contrast (right panel) for 
z = 1000 (red), 100, 10, 2, 1, and (blue). Note that the evolution of the dark energy component 
still follows the matter flow but develops an underdensity instead of an overdensity. See Section 4.4 
for more details. 



following approximations valid for \5 W \, \ a\ <C 1: 

p _ pOnt 

c 2 s ~ ato out (l - 5 W ) , 5 P = ~-a5 w . (4.9) 

From the previous relations it is easy to understand how speed of sound and pressure gradients 
are related to the density contrast. In particular it is S' p ~ c' s ~ —S' w , as one can also see 
by comparing the left panel of Fig. 9 with the right panel of Fig. 7. As a consequence, the 
pressure gradients will push the dark energy component out of the free-falling geodesic away 
from the overdense regions, as shown in the right panel of Fig. 9. One can indeed obtain at 

2 

linear order (see, e.g., [80]) the approximate relation v p ~ — 1+ ^ 8' w . 
4.4 Evolution for w < — 1 

Finally, we will now consider the case of w ou t = — 1.2 and a = — 10 -10 . We chose the same 
order of magnitude for the (positive) speed of sound of the previous Section for comparison 
sake, see the left panel of Fig. 11. Clearly, we can decrease the magnitude of a in order to 
prevent possible pathologies [62]. 

As one can see from the left panel of Fig. 10, the evolution of the matter component is 
similar to the previous case, with a moderately higher contrast at late times probably due to 
the different expansion history: for u> ou t = — 1.2 the dark-energy dominance happens indeed 
later with respect to the case for u> out = —0.8. 6 Quite different is instead the evolution of 
the dark energy component, which still follows the matter flow but develops an underdensity 
rather than an overdensity (see the right panel of Fig. 10). This can be understood from 
the fact that the relation 5 w /5m ~ /jT'™' changes sign for w out < — 1. It is interesting 
to note that a similar behavior was found in the analysis of Ref. [66] of perturbations in 
scalar-tensor cosmologies, where the non- minimal coupling of the field is indeed responsible 

6 Similar considerations apply to the case of a cosmological constant which gives a present-day contrast of 
about ~ 8, roughly half way between the cases studied in Sections 4.2 and 4.4. 
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Figure 11. Left panel: evolution of the speed of sound for the times corresponding to the 
background redshifts of z = 1000 (red), 100, 10, 2, 1, and (blue). Right panel: evolution of the 
peculiar velocities of the dark-energy component with respect to the matter component. See Section 
4.4 for more details. 



for the crossing of the so-called phantom divide (w ou t = — 1). As one can see from the right 
panel of Fig. 11, the evolution of the peculiar velocities is instead approximately the same 

C 2 I 

one of the previous Section, as in the relation v p ~ — S w both S w and 1 + u> ut change 
sign. 

5 Comparison with previous work on exact solutions 

In this Section we will compare our findings to previous work dealing with exact solutions. As 
explained in Section 1 the model here presented is a generalization to the case of n decoupled 
and non-comoving perfect fluids of the Lemaitre model [1] (see [2-5] for recent contributions) 
which describes the dynamics of a spherically-symmetric perfect fluid. Therefore, the model of 
Eqs. (2.32-2.36) straightforwardly reduces to the Lemaitre model, possibly with the addition 
of a cosmological constant. 

In the case of one source only, our equations can be used to describe the fluid in non- 
comoving possibly- geodesic coordinates, which may be desirable if the fluid has pressure 
gradients and its four-velocity undergoes acceleration. Within our formalism this is accom- 
plished by simply taking a geodesic reference- frame four- velocity it" f , i.e., by setting A' = 
(see Eq. (A. 18)). This can be useful to disentangle the effect of pressure gradients on the 
evolution of the model as shown, for instance, by Eq. (2.17) for the time evolution of the 
curvature function. In the latter equation, indeed, only the second term on the right-hand 
side contributes if comoving coordinates are used and only the first term on the right-hand 
side contributes if geodesic coordinates are used. Spherical models have been discussed pre- 
viously in the literature in non-comoving coordinates, albeit with a different purpose. The 
main idea was indeed that an exact solution with a simple appearance in non-comoving co- 
ordinates may become extremely complex when transformed to the appropriate comoving 
system, or that the integration of the partial differential equations necessary to obtain the 
comoving coordinates may even represent an intractable mathematical problem. We refer 
the interested reader to [81-87] and to [88] and references therein. 
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Another line of research focused on interpreting two or more fluid components as a single 
effective fluid, see for example [89] and references therein. The work of [90, 91] showed indeed 
that a mixture of fluids whose four- velocities lay on a two- plane can always be reinterpreted as 
a single anisotropic fluid. The model of this paper satisfies this last requirement and it would 
be interesting to investigate in future work the amount of anisotropic pressure generated by 
the evolution of a cosmological multi-component fluid modeling, e.g., baryons, dark matter 
and dark energy. As a simple example we can apply here the work of [90] to the model of 
Section 3 featuring two non-comoving dust components in a flat ACDM universe. Following 
[90], the energy- momentum tensor of the effective fluid is: 

r e f = p eff n>f ff + 5 e f , (5.1) 
and a short calculation shows that the rest energy density is: 

Peff = \(PM + PN) + ^V(PM + PN) 2 + 4p M pN (l 2 ~ 1) , (5-2) 

and that the components of the diagonalized stress tensor are: 

PA = 0, (5.3) 

PR = -l^iPM + PN) + ^V(PM + PN) 2 + 4p M PN {l 2 ~ 1) , (5-4) 

where pa and pr are the angular and radial pressures, respectively. Clearly, for 7 — > 1 
(v c — > 0) the effective fluid becomes simply the sum of the two fluids and Pa = PR = 0. For 
a nonzero radial velocity v c 7^ 0, however, the radial pressure is positive pn > and the 
effective fluid is anisotropic. See [90] for more details. Finally, we would like to mention that 
the authors of [91] claim that the opposite process is also feasible, i.e., that given a particular 
anisotropic fluid, it is always possible to find a multi-fluid model with an equivalent energy- 
momentum tensor. Note, however, that this last claim was argued not to be valid for a 
general anisotropic fluid [92] such as the anisotropic fluid models of compact objects like 
neutron stars [93, 94]. 

Finally, [95] studied anisotropic pressure in the context of the LTB metric in order 
to investigate the evolution of inhomogeneities in an interacting photon-baryon mixture. 
A physically-plausible hydrodynamical description of cosmological matter in the radiative 
era between nucleosynthesis and decoupling was given, emphasizing its thermodynamical 
consistency. It would be interesting to develop such considerations also in the presence of a 
dark matter component using the multi-fluid model presented in this paper. 

6 Conclusions 

We have extended the Lemaitre model to the case of n decoupled and non-comoving perfect 
fluids with general equations of state. We have expressed the full set of 3 + 2n exact equa- 
tions governing the dynamics of the model in a concise and transparent way thanks to the use 
of physically meaningful quantities like, for example, expansion and acceleration rates, ex- 
pansion and acceleration scalars, total derivatives along four-velocity and acceleration. This 
general solution can have many possible applications, both at early and late times. 

At late time, one can study the evolution of overdensities and underdensities in a uni- 
verse where dark energy is not the cosmological constant and, in particular, can be inho- 
mogeneous. This possibility may also induce an inhomogeneous variation of fundamental 
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constants if a coupling between the dark energy and the matter-radiation Lagrangian is al- 
lowed. Moreover, dark matter and baryons can be described as two separate fluids, the latter 
possibly featuring pressure. Pressure in general can have a non-negligible effect on the cos- 
mological models and this should be taken into consideration while interpreting cosmological 
datasets. At early times the contribution of radiation can be included, which may be relevant 
for the understanding of the evolution of the inhomogeneities. 

In this paper we focused on the formal understanding of the general solution, which 
we applied to two different setups. For the case of two non-comoving dust components in a 
flat ACDM universe we found rich and interesting physics as, for instance, the gravitational 
interaction between the two fluids which makes one component follow the collapse of the 
other, or the evolution of the peculiar velocities which are increasing with respect to the 
background but decreasing between the two fluids. 

Finally, for the case of clustering dark energy, we have been able to follow in an exact way 
the collapse of an overdensity to which also the dark energy contributes, and we confirmed 
that the evolution agrees with the theoretical growth rates during matter domination. We 
also considered the case of a small speed of sound (c s ~ 1CT 5 ) and of a phantom equation 
of state, which also showed the expected dynamics. These applications show some of the 
interesting features of this n-fluid component solution, which we will further investigate in 
forthcoming work, both from a theoretical and observational point of view. 
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A Expansion tensor 

In this Appendix we will study the kinematical properties of a generic four- velocity field u a , 
which may be identified with the r.ih fluid-component velocity uf or with the reference- frame 
velocity u" f . 

The acceleration of u a is given by a a = u^V ^u a , whose components are: 

at = —e x ^v p a and a r = 7 — a , (A.l) 

Vc 

where the acceleration scalar a is 

a = ^jg^a a a p = v^ 1 (& T + v 2 p 9r) , (A.2) 

where Qt and Or are defined below. The expansion tensor O a p is then: 

@ a = WJfp = V^u^ + a {a u p) , (A.3) 

where h a p = g a g + u a up is the projection tensor on the hypersurface orthogonal to u a . 
In the second equality we have used the fact that u a V pu a = 0. As we are considering a 
spherically-symmetric spacetime the velocity is irrotational, i.e., there is no vorticity: 

u a p = W a h v pV {ll u v] =0. (A.4) 
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h a p is often referred to as the "spatial metric" because it is the induced metric on the space 
slices, h a p = h^ah v p g^ v . Moreover, h a p is the first fundamental form of the hypersurface and 
the expansion tensor is (minus) the extrinsic curvature, K a p = —@ a /3, which is the second 
fundamental form of the hypersurface. 

The expansion scalar can be obtained without calculating explicitly the components of 
the expansion tensor by remembering that it is: 

6 = V Q u« = = d a u« + ^ , (A.5) 

J dr 

where J = yj—g = JtJrJa, g is the determinant of the metric and we have defined the 
following quantities: 

yi 

J T = e x , J R = -j= and J A = Y 2 sin (9 . (A.6) 
V 1 + 2E 

In the last expression of Eq. (A.5), the first term can be interpreted as the "newtonian" 
expansion scalar and the second as the contribution from the metric. We then define the 
temporal, radial and angular expansion scalars: 

G T = dtu* + = e" A 7 + iv c A' , (A.7) 

CLT 

Or = d r u r + = (jv c )' + j(h r + S r } , (A.8) 

G A = d e u e + d^u* + = 7 Uh a + 2S A ) , (A.9) 



dr 

respectively, where we defined the quantities: 

Y" E' Y' 

Sr = Vc Y 7 ~ Vc \ + 2E and S A = v c y, (A.10) 

which are defined analogously to Hr and H A of Eqs. (2.2-2.3) and can be understood as 
"spatial" expansion rates relevant for non-comoving fluids. From Eq. (A. 2) follows that the 
expansion scalar for a geodesic four- velocity field is geo d. = ©_r/7 2 + ©A- 

In order to calculate the shear we need the non- vanishing components of Q a p: 

v 2 

@ a = e 2 W4 (6 T + Or) , & r t = -e A 7 2 - (©t + &r) , (A.ll) 

y Vc 

2 

&rr = 1 2 \{Qt + Qr), Q ee = @ H /sm 2 e = lY 2 Q A . (A.12) 

Note that ®tt and Q r t are nonzero as we are using coordinates in general not comoving with 
u a . The components of the expansion tensor clearly satisfy: 

e = g afi e a/3 = e T + e R + e A . (A.13) 

The shear tensor a a p is then defined in terms of G a( g as 

CT a fB = ©a/3 ~ ^© Kfl , ( A - 14 ) 
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and its non-vanishing components are: 



2 n \ o n 2, \ nV n 



a tt = ^ e2 V^ <r , <?rt = --j=e x ~f 2 ^ a , (A.15) 



2 i£ a. 



where the invariant a of the shear tensor is given by: 

a 2 = Ig^g^a^a^ = Ue T + Q R - 2 . (A.17) 

We may now particularize the above expressions for the case of the reference-frame 
velocity u" { : 

a r f = — A' , a r r f = A , at r f = , (A. 18) 

rf = e^rf + o Ar{ = h r + 2h a , V3a I{ = e fiirf - ^e Arf = h r - h a . (A.19) 

Note that, as Qt vanishes, A appears in the expansion scalar and shear only through the 
derivatives with respect to the proper time d/dr r i = u^d a = e~ x d/dt. This shows that a 
non-geodesic reference frame only changes the proper time with which the expansion rate 
is measured, but not its functional form (different is the case for the acceleration rates of 
Eqs. (2.4-2.5)). Moreover, these last results show that the definitions in Eqs. (2.2-2.3) are 
relative to the reference-frame velocity. 

B The Lemaitre metric 

It is useful to consider the case of only one fluid component, which could be seen as the 
Lemaitre metric [1] in a possibly non-comoving reference frame. We will discuss some prop- 
erties of this solution, which may help us better understand the multi-fluid non-comoving 
case. We will not discuss the dynamical evolution of the Lemaitre metric, which can be easily 
deduced from the set of Eqs. (2.32-2.36). 

Let us start by considering the Eqs. (2.30-2.31) for F and F' . By taking the combinations 
giving the derivative of F with respect to r and a it is possible to obtain the following simpler 
expressions: 

where the Euclidean and actual volumes of the space slices are, respectively: 

r 47r r Y 2 Y' 

V e = 4vr / Y 2 Y'df = —Y 3 , V = 4vr / df . (B.3) 

Jo 3 Jo v 1 + 

Note that Eqs. (B.1-B.2) give the gravitating mass of a non-comoving fluid. In the fluid rest 
frame (the actual Lemaitre metric) the previous equations have the familiar form: 

F = -4irY 2 Yp = -pV e , (B.4) 
F' = AirY 2 Y' p = pVe , (B.5) 
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so that the first gives the time evolution of F, which is constant only for vanishing pressure, 
and the second links the local density to the integrated gravitating mass. Eqs. (B.4-B.5) can 
be combined into the differential expression: 

dF = p dV e | t =CO nst - P d V e | r . =CO nst , (B . 6) 

which has the usual thermodynamical interpretation for a fluid in equilibrium (at constant 
entropy) : for an increase dV of volume at constant t the total mass-energy is increased by the 
mass-energy density of the volume dV, and for an increase dV at constant r the total mass- 
energy is decreased by the work of the pressure against dV. In this picture the coordinate r 
has to be seen as a label for the mass-energy density. 

It is interesting to point out that the previous equations use the Euclidean volume 
element and so the gravitating (or Euclidean) mass F does not coincide with the invariant 
mass M which is related to the local density by [7]: 

M ' =i *7mE p=pV '- (R7) 

We can then evaluate the time derivative of M at constant r: 

Y 2 Y' AJrJapY ^_ A _ f r Y 2 Y' Jp , Xi 



M = 4vr / ; p y f 1 -f df = 4n / , p - + e A 9 rf df 

Jo VT+2E JrJaP Jo y/T+2E \P J 

= " 47r Z^ e x e ri pdf = -pV, (B.8) 

JO V 1 + It, 

where to go from the first to the second line we have used the conservation equation (2.35) 
and for the last equality we have assumed a homogeneous pressure p' = 0. Note that if p' = 0, 
then from Eqs. (2.36) and (2.33) follows that E = 0, even if M and F are evolving. Similarly 
to Eq. (B.6), we can combine Eqs. (B.7-B.8) into: 

dM = p dV | t =const ~pdV\ r =const • (B . 9) 

Note, however, that in this last equation not only a different volume is used with respect 
to Eq. (B.6), but that Eq. (B.9) is not valid for p' ^ 0, which implies a time-dependent 
curvature (that can be interpreted as total energy per unit of mass) . In this latter case the 
usual thermodynamical expression, where only the pressure at the boundary matters, does 
not seem to hold and the value of the pressure at any r matters as shown by the last integral 
in Eq. (B.8). 

Finally, it is interesting to rewrite Eq. (2.32) as: 

o 2GF 2E k. . I1Z\ 

H A = ^r + y2 = 3<p>« - { j) e » (B.10) 

where the Euclidean average for the generic quantity X is defined as: 

, v 4vr C Y 2 Y' Xdf 

(X) c = Jo .. , (B.ll) 

and 1Z is the spatial Ricci scalar of the rest frame, which is the trace of the Ricci tensor of 
the metric gij on the hypersurface of constant t: 

MFYY 

n = ~\wr- (B.12) 
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The factor of 1/6 in Eq. (B.10) can be understood by the fact that for the FLRW model it 
is 1Z/6 = k/a 2 . From Eq. (B.10) follows that the expansion rate is sourced by the Euclidean 
average of the local density and curvature, and not by the actual averages. This is a potential 
source for "strong" backreaction effects [96, 97] of the inhomogeneities on the evolution of 
the background (see, for example, Appendix B of [40] and also [98, 99]). 

C The W n and 5 n functions 

The W n (x,a) step functions introduced in [48] are C n everywhere and go from 1 to for x 
from to 1 with a transition which is steeper for higher a G [0, 1[. In this paper we used W\ 
and W3, which have the following explicit form: 



Wi(x, a) 



\ + \ sin 



7T I \ 



x—a 
1-a 



< X < a 

a < x < 1 

1 < X 



and 



W 3 (x,a) 



1 

4tt 2 
1 

471- 2 





1+7T 2 



4-8[f=2 

' 1 — ot 



-l + 8vr 2 (f5f 



COS (471-^ 



cos ( 4tt^—^ 



< x < a 



a < x < 



l+a 



±& < x < 1 



1 < X 



From the W n function we then built the following general contrast 5 n : 



S n (r,r t ,A l! 2,ai ! 2 



A 2 + A!^ n ( i, ai ) 0<r <r t 



A 2 W n 




ft < r < r b 

n<f 



(C.l) 



where A = Ai + A2 is the contrast between the center (r = 0) of the spherical inhomo- 
geneity and the border (r = r&), and A2 is the contrast at r t . We have used 5 n to model 
inhomogeneities in pm, -Fjv and v c in Section 3. 
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